%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code generates the Figures 5,6 & 7 for "Self-Fulfilling Prophecies, 
% Quasi-Non-Ergodicity & Wealth Inequality by
% Jean-Philippe Bouchaud and Roger E. A. Farmer This version by R.E.A.
% Farmer:  August 2 2022
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
rng(2);                       % Set rng(2) to replicate the figures in the
% paper
fnum = 1;                     % First figure number. Change this to compare 
                              % figures across diferent values of alpha
alpha = 1;                    % Change this to check for robustness
                              % Default is 1


LB = 7; UB = 12;              % Set bounds for tail estimate

N =1000000;                   % This is the number of people
%N = 100000;
per =52;                      % This is the number of subdivisions
                              % in a year
T = 300*per;                  % This is the number of periods
n = min(250,T);               % This is number of observations

gap = floor(T/(n-1));
nums = floor(linspace(1,T,n));  % These are the points where we record
% the wealth distribution
bet = 0.97.^(1/per);            % bet is the annual discount rate
epsw = 1;                       % Normalize the endowment to 1
d=50;                           % This is expected life in years
del =1- exp(-1/per/d);

lam = sqrt(del/alpha);          % Set delta and lambda


Burn = floor(T/5);              % Set the burn-in period

% Initialize
Herf  = NaN(1,n);
LW    = NaN(N,n);

z = rand(N,1);
P(:,1)  = z(:,1);
D(:,1)  = zeros(size(P));
Pimp    = NaN(1,T);
Pact    = NaN(1,T);
state   = NaN(1,T);
LQ      = NaN(1,T);
Pimp(1,1) = mean(P(:,1));
Pact(1,1) = Pimp(1,1);
state(1,1) = 1;

LQ(1,1) = log(bet) + log(Pimp(1,1));
H = epsw/(1-bet*(1-del));
LW(:,1) = ones(N,1)*log(H);

P_temp  = P(:,1);
LW_temp = LW(:,1);

tic;
for t = 1:Burn-1
    x = rand(N,1)>del;
    nx = sum(x);
    z = rand(N,1);
    s = rand<mean(P_temp);
    P_temp = x.*((1-lam)*P_temp + s*lam) + (1-x).*z;
    LQ_temp = log(bet) + log(sum(P_temp.*x.*exp(LW_temp))) ...
        - log(N) - log(H);
    LW_temp = real(x.*s.*(log(P_temp) + log(bet) + LW_temp - LQ_temp) +  ...
        x.*(1-s).*(log((1-P_temp)) + log(bet) + ...
        LW_temp- log(bet-exp(LQ_temp)))...
        + log(H).*(1-x));
end

P(:,1) = P_temp;
LQ(1,1) = LQ_temp;
Pact(1,1) = mean(P_temp);
LW(:,1) = LW_temp;
BW = exp(LW_temp);
Herf(1,1) = sum(BW.^2)./(sum(BW).^2);
st = toc;
disp(['Burn in time ' num2str(st) ' seconds']);

tic;
next =1;
disp(['Iteration ' num2str(1)]);
while next  < n
    for t = 1:T-1
        x = rand(N,1)>del;
        nx = sum(x);
        z = rand(N,1);
        s = rand<Pact(1,t);
        state(1,t) = s;
        P_temp = x.*((1-lam)*P_temp + s*lam) + (1-x).*z;
        Pact(1,t+1) = sum(P_temp)/N;
        LQ(1,t+1) = log(bet) + log(sum(P_temp.*x.*exp(LW_temp))) ...
            - log(N) - log(H);
        Pimp(1,t+1) = exp(LQ(1,t+1))/bet;
        LW_temp = real(x.*s.*(log(P_temp) + log(bet) + LW_temp - ...
            LQ(1,t+1)) +  x.*(1-s).*(log((1-P_temp)) + log(bet) ...
            + LW_temp- log(bet-exp(LQ(1,t+1))))+ log(H).*(1-x));
        if t == nums(next)
            next = next+1;
            LW(:,next) = LW_temp;
            BW = exp(LW_temp);
            Herf(1,next) = sum(BW.^2)./(sum(BW).^2);
            disp(['Iteration ' num2str(next)]);
        end
        
    end
end
%%
st = toc;
disp(['Run time ' num2str(st) ' seconds']);
BW = mean(exp(LW),2);
LBW = log(BW);
LWq= quantile(LBW,[0.5 0.9 0.99]);
alph = del/lam^2;
u = linspace(0,1,1000)-1/2;
y = gamma(2*alph)/gamma(alph)^2*(1/4 - u.^2).^(alph-1);
y(y==inf)=NaN;


%  Estimate the CDF & compute tail parameter estimate
[F,X] = ecdf(LBW);                  
Y = log(1-F);
Y = Y(isfinite(Y));
X =  X(isfinite(Y));

X1 = X(X>=LB & X<= UB);
M = size(X1,1);
Y1 = Y(X>=LB & X<=UB);
Z = [ones(M,1) X1];             % Used to estimate tail coefficient
B=Z\Y1;
tail = B(2);
disp(['Tail Slope = ' num2str(tail) newline]);

%
%v Compute some statistics

q = 100*quantile(exp(LBW),[0.5,0.99, 0.999])/(H);
Hnorm = (mean(Herf)-1/N)/(1-1/N);
Hm = mean(Herf);
%     Vm = mean(StdD.^2);
F1 = ((1:N)/N);
Ws = cumsum(sort(BW,'ascend'));
Lorenz = Ws/Ws(end);
AreaB = sum(Lorenz/N);
G = 1-2*AreaB(end);

disp(['Wealth Distribution' newline '50''th percentile = ' num2str(round(q(1),0)) '% of Human Wealth' newline ...
    '99''th percentile = ' num2str(round(q(2),0)) '% of Human Wealth' ...
    newline '99.9''th percentile = ' num2str(round(q(3))) '% of Human Wealth'...
    newline 'Human Wealth = ' num2str(round(H,0)) newline ...
    'Normalized Herfindhal Index = ' num2str(round(10^3*mean(Hnorm),4)) ' x 10^-3' newline ...
    'Gini Coefficient ' num2str(round(G,2))]);
%%
% Draw Graphs

Pls =   1:T;
% Construct PE ratio
pE = bet/2*((2*Pimp(Pls)-1)./(1-bet*(1-del)) + 1/(1-bet))/52;
% Construct Risky Return

%%

RR = ((pE(2:end-1)+(state(2:end-1)./52))./pE(1:end-2));
Rs = ones(size(RR))./(bet);
J = floor(T/52);
RR1 = NaN(1,J);
Rs1 = NaN(1,J);
for j = 1:J-1
    LRR = log(RR);
    RR1(j) = mean(52*LRR( 52*(j-1)+1:52*(j-1)+52 ));
    Rs1(j) = -52*log(bet);
end
RR1 = RR1(1:end-1);
Rs1 = Rs1(1:end-1);

figure(fnum);
subplot(2,2,1);
plot(u+1/2,y,'LineWidth',4);
title('Stationary Measure','FontSize',16);
xlabel('Value','FontSize',16);
ylabel('Probability','FontSize',16);
ylim([0,1.1]);

subplot(2,2,2);
plot(Pls/per,100*(Pact(Pls)-Pimp(Pls)),'LineWidth',2);
title('Actual Minus Implied Probability','FontSize',16);
xlabel([num2str(ceil(T/per)) ' Years of Weekly Data'],'FontSize',16);
ylabel('Percentage Difference','FontSize',16);

subplot(2,2,3);
histogram(LBW,100,'BinLimits',[-5,12]);
title({'Log Wealth Distribution'},'FontSize',16);
xlabel('Log of Wealth','FontSize',16);
ylabel('Count','FontSize',16);

subplot(2,2,4);
plot(Pls/per,pE(Pls),'LineWidth',2);
title('Price Earnings Ratio','FontSize',16);
xlabel([num2str(ceil(T/per)) ' Years of Weekly Data'],'FontSize',16);
ylabel('Price of Equity','FontSize',16);



% plot(Pls/per,exp(LQ(Pls)),'LineWidth',2);
% title('Price of a Security in the High State','FontSize',16);
% xlabel([num2str(ceil(T/per)) ' Years of Weekly Data'],'FontSize',16);
% ylabel('Security Price','FontSize',16);

d = -1/log(1-del)/per; Lm = -1/log(1-lam)/per;
sgtitle([{['Expected Life = ' num2str(round(d,2))...
    ' years']}; {['Expected Belief Decay Rate = ' num2str(round(per*Lm,0)) ...
    ' weeks']}],'FontSize',16);

%%
fnum = fnum+1;
figure(fnum);
plot(F1,[100*Lorenz'; 100*F1]','LineWidth',2);
title('Lorenz Curve for Simulated Data','FontSize',16);
xlabel(['Percentile ' newline ...
    'Gini Coefficient = ' num2str(round(G,2))] ,'Fontsize',16);
ylabel('Cumulative Wealth ','Fontsize',16);
%%
fnum = fnum+1;
figure(fnum); 
Left_plot = 0;
X2 = X(X>Left_plot);
Y2 = Y(X>Left_plot);
numpts = ceil(numel(X2));
inter = ceil(numel(X2)/numpts);
ind = 1:inter:numpts;


X_plot = X(X>LB*0.9);   % Used to plot tail coefficient
M0 = size(X2,1);
Z0 = [ones(M0,1) X2];


scatter(X2(ind),Y2(ind),'LineWidth',0.01);     % Plot graph to estimate tail
hold on;
line([LB,LB],[min(Y2),max(Y2)],'Color','black','LineStyle','--'...
    ,'LineWidth',2);
line([UB,UB],[min(Y2),max(Y2)],'Color','black','LineStyle','--'...
    ,'LineWidth',2);
plot(X2,Z0*B,'LineWidth',4);
hold off;
ylim([min(Y2),max(Y2)]);
title(['Tail Slope Coefficient = ' num2str(round(tail,1))],'FontSize',16);
ylabel('Log(1-CDF)','FontSize',16);
xlabel('Log Wealth','FontSize',16);
%%

fnum = fnum+1;
plot([Rs1;RR1]','LineWidth',3);
[mean(Rs1);mean(RR1)]





